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Abstract: Urbanization research is of great significance for the sustainable use of regional land resources, ecological 
environment protection. In the middle reaches of the Yellow River (MRYR), the expansion process and driving factors 


of urban construction land at different scales are not yet clear. This study explored the distribution pattern of 
urban construction land on different slope gradients at different scales and analyzing its influencing 
factors. The main findings were as follows: (1) Significant expansion of urban construction land 
occurred in the MRYR over the past 20 years. Spatial heterogeneity was observed in the regional 
urban construction land expansion process among different geomorphic regions. (2) The urban 
construction land in the MRYR expanded vertically to areas with slopes »5?, especially in 2005-2010. Significant 
slope climbing of urban construction land were observed in the loess hilly-gully and rocky 
mountain areas. (3) 68.45% of the counties of the MRYR belonged to the slope-climbing types, 
included 37.38% belonged to the high-slope-climbing types. (4) The regional population density and 
economic development level were closely associated with regional urban construction land area 
variability. (5) The climbing process of regional urban construction can effectively alleviate 


farmland encroachment and pressure on the regional ecological environment. The urban expansion of 
the metropolitan distribution areas in the Plain region (such as Xi 'an, Taiyuan, etc.) had a more significant impact on 
the local carbon storage. 
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1. Introduction 
Urbanization affects the evolution of land spatial patterns and functions, and has a significant impact on regional 


ecosystem stability, food security, carbon emission, and so on (Hu et al. 2023). In recent years, China's rapid 
urbanization process has intensified the supply-demand contradiction of available land resources (Zhou et al. 2022; 
Zhang et al. 2022; Liu et al. 2023), causing changes in the surface water and thermal environment, bringing enormous 
pressure to regional ecosystems (Mohammad et al. 2021; Song et al. 2022; Shi et al.,2023), and also posing a certain 


threat to regional food security (Mohammad et al. 2021). Besides, Land use change brought by urbanization 
development caused a large loss of carbon storage in the terrestrial ecosystem, which has a significant impact on the 


carbon cycle of the terrestrial ecosystem, especially in ecologically fragile regions (Gong et al. 2023; Wen et 
al. 2024; Yuan et al. 2024). A systematic and in-depth analysis of the urbanization process and the development and 
evolution of urban construction land is crucial for achieving intensive use of land resources and high-quality regional 
development. It is of great significance for sustainable regional economic development, ecosystem stability, food 


security, and achieve the dual carbon target (Zhang X et al. 2023; Zhang Y et al. 2023; Wen et al. 2024). 

The expansion of urban land mainly takes two forms: horizontal expansion and vertical growth (Zhou et al. 2021a). 
Previous researchers have conducted extensive research on the spatial and temporal distribution pattern of 
construction land (Zhang et al. 2020; Zhang Y et al. 2023), landscape evolution (Yang D et al. 2022; Gilman et al. 
2023), driving factors (Ouyang et al. 2021; Monsted et al. 2023; Mou et al. 2023), simulation prediction (Arifin et al. 
2023; Liao et al. 2023), and other aspects of regional urban horizontal expansion; And it is pointed out that there is a 
changing pattern of extension enhancement deceleration enhancement stability and unchanged" in the horizontal 
expansion of urban areas (Zhan et al. 2020), which is mainly influenced by factors such as urban development planning 
and policies, regional population size, economic development level, natural environment, transportation, etc. (Ouyang 
et al. 2021; Mou et al. 2023). Some studies have also pointed out that the rapid urban expansion leads to landscape 
fragmentation and spatial differences in agglomeration, which exacerbates the complexity of regional landscape 
morphology (Yang D et al. 2022), changes the regional water and thermal environment, and puts certain pressure on 
regional ecosystem stability and food security (Mohammad et al. 2021; Song et al. 2022). On this basis, some 
researchers have simulated the future urban expansion process using methods such as Cellular Automata and Markov 
chains, and proposed corresponding countermeasures based on the simulation results (Arifin et al. 2023; Liao et al. 


2023). The above research has greatly promoted regional urbanization research and provided guidance and experience 
for the sustainable development of regional urbanization. 

In recent years, construction land has gradually expanded to relatively higher-slope regions in some parts of the world 
owing to the shortage of land resources available for urban construction (Zhou et al. 2021a). This is particularly 
significant in complex terrain areas with large slope changes, and there has been little research on urban slope climbing 
in previous studies (Ouyang et al. 2021; Zhang Y et al. 2023). As an important form of vertical urban expansion, urban 
climbing has became a widely concerned issue in the process of urbanization in non plain areas (Zhou et al. 2021a). The 
slope spectrum, as an important indicator for quantitatively studying the process of urban climbing, can quantitatively 
describe the combined morphological characteristics of slope surface in terrain gradients through spectral lines, 
including terrain slope spectrum curve and construction land slope spectrum curve (Wang et al. 2007; Zhou et al. 
2021b). In recent years, some scholars have conducted research on urban slope climbing processes under different land 
use modes using slope spectra (Zhou et al. 2021b; Jiang et al. 2022; Zhang H et al. 2023). By constructing the Average 
Built-up Land Climbing Index (ABCD), researchers systematically analyzed the changing characteristics and patterns of 
construction land slope spectra in China during 1990-2018 and pointed out that the climbing phenomenon is significant 
in the complex terrain of the central and western regions, and there is significant provincial heterogeneity in the degree 
of construction land climbing (Zhou et al. 2021b). Besides, some researchers explored the climbing process of 
construction land in typical cities such as Shenzhen (Peng et al. 2018), Guiyang (Peng et al. 2023), and Lanzhou-Xining 
(Zhang H et al. 2023). They found that the phenomenon of urban construction land climbing occurred in all of these 
cities. Among them, the scale and degree of urban construction land climbing in the Lanzhou-Xining urban cluster 
showed a gradually decreasing trend (Zhang H et al. 2023). However, current research on urban expansion based on 
slope range mostly focuses on urban research under a single landform type (such as valley towns (Jiang et al. 2022), 
mountain towns (Peng et al. 2022), low hills and gentle slope settlements (Peng et al. 2015), etc.), lacking systematic 
and in-depth research on vertical urban expansion under complex landform environments (Zhou et al. 2021b; Zhang H 
et al. 2023), especially in the middle reaches of the Yellow River (MRYR) where the ecological environment is fragile 
and the contradiction between humans and land is sharp, the process and driving factors of the expansion of urban 
construction land in this region are not yet clear. The influence of urban expansion on the regional environment and 
carbon reserves still needs to be further explored. Furthermore, previous studies on vertical expansion of urban 
construction land mainly focus on the scale of country (Zhou et al. 2021b), urban agglomeration (Zhang H et al. 2023) 
and city (Peng et al. 2018), and rarely involves the study on the county scale. However, there is a significant spatial 
heterogeneity in the level of urbanization in Chinese counties, especially in the complex terrain of the MRYR (Dai et al. 
2024). Therefore, it is necessary to carry out relevant research on the vertical expansion of urban construction land 
from the county level. 

The MRYR is located in a transitional zone between arid and semi-arid regions in China, with complex terrain, an arid 
climate, and prominent ecological and environmental problems. The rich and diverse mineral resources in this region 
make it the primary energy production base for coal and other energy sources in China (Zhao et al. 2021). Recently, 
with the large-scale development of coal resources and rapid growth of the regional population, the level of urbanization 


in the MRYR has significantly improved and the contradiction between supply and demand of land 
resources increased (Fan et al. 2022; Wang et al. 2022a; Qiao et al. 2023). 

Therefore, this study focuses on urban construction land in different landforms in the MRYR and systematic 
analysis of the spatiotemporal distribution pattern of urban construction land on slope gradient at 
the geomorphic and county scales. The degree of influence of comprehensive regional 
environmental factors on urban expansion in different terrains was explored using a structural 
equation model. The influence of urban expansion on farmland, ecological land and regional carbon 
reserves was also evaluated. The results provide a theoretical reference for regional urban 
construction land planning, sustainable utilization of land resources, and alleviation of the 
contradiction between humans and land. 


2. Data Sources and Methodology 

2.1 Study area 

The MRYR is located at approximately 33°45'-40°11'N and 104?27' -113?39'E, adjacent to the Taihang Mountains to the 
east, Liupan Mountains to the west, Qinling Mountains to the south, and Yinshan Mountains to the north (Jia et al. 
2023). It has a length of 1,206.4 km and a basin area of 344,000 kme, accounting for 43.3% of the Yellow River 
basin area, and mainly includes the six provinces and autonomous regions of Shaanxi, Shanxi, Henan, Gansu, Inner 
Mongolia, and Ningxia (Fig. 1a). Most areas have a temperate monsoon climate, which is hot and rainy in summer and 
cold and dry in winter. Precipitation is concentrated from June to September and is unevenly distributed across 
different regions (Gou et al. 2019). From 2000 to 2020, the average annual precipitation in the region was 
between 400 and 600 mm, and the average annual temperature was between 9 and 10'C (Fig. 1b). The 
MRYR is one of the most vulnerable areas in China because of its sparse vegetation, apparent terrain cutting, and 
interaction between wind and water erosion. This region has a long history of agricultural civilization, 


with a relatively weak economic foundation and concentrated population distribution. In recent 
years, the level of regional urbanization has led to rapid urbanization due to rapid population 
growth brought about by the exploitation of mineral resources. By 2021, the population of the 
MRYR was more than 19.8 million, accounting for 13.96% of the total population of China, and 


the gross national product was projected as 13179.28 billion yuan, accounting for 11.52% of the national 
GDP (Guo et al. 2024). 

The MRYR is characterized by complex topography, which includes ravines. To further analyze the urbanization process 
in different geomorphic regions of the region, the research area was divided into five categories according to natural 
conditions (Jin et al. 2012): the northwest sandstorm (SA), loess hilly-gully (HGA), loess plateau-hilly (PHA), plain 
(PA), and rocky mountain regions (RMA). 

2.2 Data sources 

The data in this study mainly included remote sensing, land use, and meteorological data (Table 1). This study used 
impermeable surface data to reflect the changes in urban construction land in the study area (Qian et al. 2019). Land 


use and cover maps of the MRYR for five periods (2000, 2005, 2010, 2015, and 2020) were obtained from the 
Annual China Land Cover Dataset (CLCD, https://zenodo.org/records/5816591) with a resolution of 30mx30m, and 
the primary land categories included cropland, forest, grassland, shrubland, water, snow/ice, impervious, 
barren, and wetland. According to previous studies, the land use data of the CLCD has relatively 
higher spatial resolution and classification, which meets the requirement of the research precision 
(Yang et al. 20212; Liu et al. 2023). Nighttime lighting data were used to reflect regional economic 
activities and per-capita income (Gibson et al. 2021). The nighttime lighting data processing and principles are 


referenced from previous studies (Chang et al. 2023; Zheng et al. 2024). The details description was shown in 
Appendix 1. According to the requirements of the Code for Vertical Planning of Urban and Rural Construction Land 


(CJJ83-2016) and regional topographic slope spectrum of the MRYR, and considering the suitability of the 


construction land, this study only analyzed construction land with a slope «25? (Zhang etal. 2023). The general 
framework of this study was shown in Fig.2. 

2.3 Method 

2.3.1. Urban expansion spatial distribution research methods 

In this study, a spatial analysis of urban expansion was conducted using the topographic slope spectrum curves and 
construction land slope spectrum curves (Wang and Tang, 2007; Zhou et al. 2021b). The slope data were rasterized with 
the 5-period land-use data (2000,2005,2010,2015 and 2020), the number of rasters for each slope interval and the 
number of regional construction land rasters were counted at 1°intervals, and then the topographic slope spectrum 
curves and construction land slope spectrum curves were drawn. 

The topographic slope spectrum (land area frequency slope range curves) is a statistical map of the percentage of the 
land area corresponding to the slope factor in a specific statistical area to the total statistical area, which can visually 
reflect the overall slope distribution of the background land (Zhang et al. 2023). The construction land slope spectrum 
(construction land area frequency slope spectrum) is a statistical map of the percentage of the construction land area 
corresponding to the slope factor in a specific statistical area to the total construction land area in the statistical area, 
which can visually reflect the overall slope distribution of the background land and the land-use preference of 
construction land at different slopes (Zhang et al. 2023). 

a. The topographic slope spectrum and the construction land slope spectrum 

The topographic slope spectrum (i.e., land area frequency slope spectrum) and the construction land slope spectrum 
(Zhou et al. 2021b; Zhang et al. 2023) were calculated as follows: 
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where P; and Pci, are the topographic slope range and construction land slope range of x raster cells, respectively; 
Area:x and Aredcix are the land area and construction land area of x raster cells, respectively; and Arearand Areaci are 
the total land area and construction land area of the study area, respectively. 
When the Perx=Ptx, the slope range of construction land and the ground slope range intersect, and the slope at this time 
is defined as the advantageous slope threshold (T). When Pc» P, the slope range of construction land with slopes less 
than T is larger than the terrain slope range, and construction land is mainly distributed in the raster areas with slopes 


less than T at this time (called as the advantageous area). Conversely, the other cases are defined as the inferior 
area (Pax « Pi). 


b. Average built-up land climbing index and upper limited slope 
In this study, the average built-up land climbing index (ABCT) (Zhou et al. 2021b; Zhang et al. 2023) was used to 


quantitatively measure the degree of construction land climbing. The specific calculation formula 
is as follows: 


BCI= Sa — a )x100% (3) 
ABCI = BCI/(b-a) (4) 


where BCI (defined as the proportion of the construction land area in the disadvantaged area to the total construction 


land area in a given period of time) is the built-up land climbing index from year a to year b; GAa and GA; 
represent the areas of construction land in the disadvantage area in year a and year b, respectively; Aa and Ap are the 
total areas of construction land in year a and year b, respectively; ABCI is the average climbing index of construction 


land from year a to year b. ABCI>0 indicates that the proportion of construction land in the disadvantaged area 
increases with time (construction land exhibits climbing expansion). The larger the ABCI value 
respects the stronger the climbing capacity of the construction land. 

The upper limited slope (ULSC) (Zhou et al. 2021b; Zhang et al. 2023) was calculated to analyze the difference 
in slope maxima over time for different construction land distribution scales. When ULSC > 0, it 
indicates that the construction land is climbing to the disadvantaged area in that period, and the 
larger value of the ULSC means the greater the climbing capacity of the construction land. ULSC 
< 0 indicates that the construction land is expanding horizontally in the advantaged area in that 


time. 


2.3.2 Regulation effect of terrain on urban construction land research method 

Structural equation modeling (SEM) is a multivariate statistical method based on a covariance matrix that includes 
factor and path analyses and other statistical methods (Gu et al. 2022). It can be used to analyze the relationship 
between latent and observed variables, as well as the cause-effect relationship in the latent variables. The path 


coefficient beside the arrow in the plot drawing by the SEM respect the degree of influence between variables or 
the degree of influence of variables (Wang et al. 2021). The path coefficients range from -1 to 1, with values closer to 1 or 
-1 indicating a stronger relationship between the variables and a coefficient of o indicating no relationship. The total 
effect of one variable on another is equal to the sum of the direct and indirect effects. The direct effect is the path 
coefficient of the variable pointing directly to the target variable, whereas the indirect effect is the product of the path 
coefficients pointing from the variable to the target variable (Xiong et al. 2023). 


In this study, the AMOS (Analysis of Moment Structure) in SPSS software was used for SEM analysis, 


and the Pearson correlation analysis was used to analyze the significance. Based on SEM, this study 
explored the impact mechanisms of multiple factors on urban construction land, including temperature (Tem), 
precipitation (Pre), population density (Pop), nighttime lighting (NTL), normalized vegetation index (NDVI), and 
potential evapotranspiration (PET), within a specific range of slope, slope aspect, and altitude. 


2.3.3 Carbon storage calculation model 
The carbon storage module of the INVEST model was used to analyze the carbon storage variation in the MRYR 
between 2000 and 2020, which classified the carbon density of each land use type as four primary carbon pools: 
aboveground living biomass (Cabove), belowground living biomass (Chelow), soil carbon (Csoi), dead organic matter (Caeaa). 
According to the distribution of land use types in the MRYR, the InVEST model analyzes and calculates the total carbon 
storage (Ctotal) of various land use types (Shi et al. 2024; Wang et al. 2023; Li et al. 2023). The specific formulas are as 
follows: 
Ctotal=Cabove+ Cbelow-- Csoil-- Cdead (5) 

Based on carbon density and land use data from different regions, the formula for calculating carbon storage for each 
land use type is: 

Ctotal= (Cabovet+Cbelow+ Csoil+Caead) x Ai (6) 
i represents the average carbon density of each land use type, and Ai represents the area of that land use. The total 
carbon storage is the sum of all carbon storage. 
In this study, the classification and carbon pool of each land use type were determined by summarizing previous 
research results (Wang et al. 2022b; Yang et al. 2021b; Duan et al. 2024), as shown in Table 2. Parameters of Carbon 
density in this study has been corrected (annual average temperature and biomass carbon density were corrected by the 
formula from Giardina et al (2000) and Chen et al (2007), regional precipitation factor was corrected from Alam et al. 
(2013); Detailed formulas are described in Appendix 1). 


3. Results and Analysis 

3.1. Overall characteristics of the environment, population, and economy in the MRYR since 2000 

During the period 2000-2020, relatively smaller fluctuations of the annual mean temperature were observed in the 
MRYR, and a gradually decreasing trend of annual mean temperature appeared from southeast to northwest (Fig. 3a). 
Higher air temperatures were observed in the PA and RMA. Annual precipitation in the MRYR showed a gradually 


increasing trend and displayed gradually decreased from southeast to northwest in the MRYR (Fig. 3b). Regional PET 


presented a spatial distribution pattern of gradually decreasing from southeast to northwest. From 2000 to 2020, 


relatively higher NDVI values were observed in PA and RMA. Overall, vegetation coverage in the study area increased 
significantly, especially during the 2010-2020 period (Fig. 3d). 
During the study period, the areas of cropland and barren land in the MRYR showed a decreasing trend, while the areas 


of forest, water, and impervious surface showed an increasing trend (Fig. 3e). However, the regional grassland area did 
not fluctuate significantly during the study period. 

Overall, the regional population distribution was dense in the southeastern region and sparse in the northwestern 
region, and the population density has gradually increased in recent years (Fig. 3f). During the study period, the total 


GDP of the MRYR showed an increasing trend (Fig. 3g). Areas with higher GDP were mainly concentrated in the 
southern and eastern parts of the MRYR, particularly in the PA and RMA regions. 


3.2. Urban construction land expansion process in the MRYR 

3.2.1. Time change of urban construction land expansion in the MRYR 

The urban construction land area in the MRYR has shown a remarkable expansion trend over the past two decades (Fig. 
4a—d), gradually increasing from 7.3x103 km? in 2000 to 12.9x103 km? in 2020. The urban construction land expansion 
area was mainly concentrated in the larger cities of the Fenhe and Weihe River basins. As for the growth rate, the urban 
construction land expansion rate in the MRYR experienced a *declining-rising-declining" trend during 2000-2020, 
with a declining trend from 2000 to 2005 and an increasing trend from 2005 to 2010. A lower regional urban 
construction land area expansion rate was observed after 2010, and then a relatively stable expansion rate occurred in 
the 2015-2020 period (Fig. 4d). 

During the study period, urban construction land in each geomorphic zone showed a significant increasing trend (Fig. 
5), with considerable urban expansion occurring in the SA and HGA regions. With regional coal resource development 
and population growth, the urbanization level in SA significantly improved, and the area of urban construction land 
increased significantly (from 0.28x103 km? in 2000 increasing to 1.18x103 km? in 2020, with growth rate is 321%) (Fig. 
4). The construction land area in the HGA region has doubled in the past 20 years (9.2x103 km? in 2000 to 19.1x103 
km? in 2020), with the fastest growth rate occurring during the 2005-2010 period. The expanded urban construction 
land area was concentrated in the western and northwestern parts of the region (Fig. 5). Owing to the flat terrain and 
concentrated population distribution, most of the urban construction land in the MRYR is mainly distributed in the PA 
region. Until 2020, the area of urban construction land has expanded to 67.1x103 km? (39.1x103 km? in 2000), with an 
overall growth rate of 71.6196. The rate of increase in urban construction land area in the PA region during 2015-2020 
displayed a slowing trend. Urban construction land expansion in the PA and RMA was mainly concentrated in the areas 
surrounding the original towns. 

3.2.2. Vertical expansion of urban construction land in the middle reaches of the Yellow River 

The urban construction land in the study area presented a regular distribution in slope gradient from 2000 to 2020 
(Fig. 6a), with a skewed distribution and spectral peaks in the slope range of 0°-2°. Above a 2° slope, the proportion of 
construction land area gradually decreased and eventually tended to o. The slope spectrum of construction land 
intersected with the terrain slope spectrum at a slope of approximately 5?, indicating that regional urban construction 
land was mainly distributed in regions with slopes below 5?. The urban construction land exhibited a trend of vertical 
expansion to the areas above slopes of 5?, with the area increasing by approximately 9796 over the past 20 years.The 
urban construction land area in areas with slopes below 5? in the MRYR decreased from 92.896 in 2000 to 91.0996 in 
2020. 

The variation in the ULS of construction land in each geomorphic area also presented spatiotemporal heterogeneity, 
with obvious fluctuations in the ULS value observed in the SA regions. The urban construction land ULS value in the SA 
region showed an apparent downward trend (from 11? in 2000 to 6? in 2020) from 2000-2020, implying that the urban 
construction land in this region gradually expanded to the lower slope region. Comparing the ULSC of urban 
construction land in different periods (Fig. 6g), a relatively higher upper limit slope value (above a 1? slope) of urban 
construction land distribution was observed in the 2000-2005 period and varied weakly for the remaining periods. 
Besides, a significantly rising trend of ABCI was observed in the 2000-2010 period, followed by a sharp downward 
trend in 2010-2015. 

Based on the results of the regional slope spectrum and upper slope analysis, the peak area of urban construction land 
slope-climbing and vertical expansion mainly occurred in 2005-2010 in the MRYR and was higher than the average 
level in China during the same period (Zhou et al. 2021b; Zhang H et al. 2023). Unlike other regions in China, the 
vertical urban construction land expansion rate in the MRYR gradually trended downward between 2010 and 2015, 
which was closely related to the implementation of regional soil and water conservation and vegetation restoration 
projects after 2010 (Wang et al. 2024). 

Owing to the differences in the geographical environment and urban-rural development stages (Zhou et al. 2021b), the 
distribution and climbing degree of urban construction land in the five geomorphic areas in the MRYR differed (Fig. 
6b-f). During the study period, the peak value of the construction land slope spectrum in all geomorphic areas showed 
a downward trend, among which the peak decline in the SA region was the largest (4.2196), followed by those in the 
RMA (2.3296), PHA (0.6196), PA (0.3996), and HGA regions (0.3196) (Fig.6). The urban construction land slope spectra 
in each geomorphic type area presented a smooth inverted "V" curve. Among them, the urban construction land slope 
spectrum in the SA region was similar to regional topographic slope spectra (Fig.6). The peak value of the urban 


construction land slope spectrum in each geomorphic type area of the study area occurred on a 1° slope in 2000-2020, 
indicating that the urban construction land area in the MRYR was mainly distributed in the 1° slope region. The urban 
construction land slope spectrum in each geomorphic zone presented significant fluctuations in the 0°-5° slope region, 
while the slope spectrum curve showed a decreasing trend with increasing slope in the region of 5° to 25° (Fig. 6). The 
proportion of construction land area located below 5° slopes in all geomorphic zones exceeded the average value in 
China, except for that in the SA region (in 2000) (Zhou et al. 2021b). 


Fluctuations in the ABCI were observed in each geomorphic division during the study period (Fig. 6h). The ABCI in the 
PHA, PA, and RMA regions showed similar variation characteristics and remained stable. In contrast to other regions, 
the ABCI of SA displayed a significant upward trend during the study period, with an ABCI value «o, indicating that the 
slope-climbing of urban construction land in this region was not evident. Over the past 20 years, ABCI variability in the 
HGA region demonstrated stronger urban construction land slope-climbing, especially in the 2000-2005 and 2015- 
2020 periods. 

3.2.3. Vertical expansion of county-level urban construction land in the MRYR 

County, as an important carrier of new urbanization, is the most direct operational layer for regional land use and 
management (Sun et al. 2017; Ren et al. 2023). The overall level of urbanization in Chinese counties is low and regional 
differences are significant, especially in the complex terrain of the middle reaches of the Yellow River (Dai et al. 2024). 
Based on the ABCI and ULSC calculation results (Fig. 7), 206 counties in the MRYR were divided into three categories: 
high-climbing (ABCI > o and ULSC > 0), low-climbing (ABCI > o and ULSC > 0, or ABCI < o and ULSC > 0), and 
horizontal expansion types (ABCI < o and ULSC x o). The number of climbing counties in the study area (141) 
accounted for 68.45% of the total number of counties in the MRYR (37.38 and 31.07% high- and low-climbing types, 
respectively) during 2000-2020. The largest proportion of hill-climbing counties was in the RMA region (89.1996), 
whereas the lowest proportion was in the PHA region (37.0496). Typical high-climbing counties included Jincheng in 
the RMA region, Lingshi in the PA region, and Loufan County in the HGA region, among which Jincheng had the 
highest ABCI (1.7496). The low-climbing counties included Baode in the HGA region, Wenxi in the PA region, and Uxin 
Banner in the SA region (Fig. 7). The horizontally expanding counties were mainly concentrated in Xi in the HGA 
region, Jingchuan in the PHA region, and Heshun in the RMA region. 

The increasing degree of county urban construction land in the MRYR showed significant spatial-temporal 
heterogeneity (Fig. 7). Owing to regional economic development and the exploitation of mineral resources (Wu et al. 
2022), higher urban construction land slope-climbing was evident in the SA region from 2000 to 2005.Although flat 
land resources in the PA region were relatively abundant, construction land in this region gradually expanded toward 
areas with higher-slope degrees due to the need for regional economic development and population survival (Peng et al. 
2023). In the RMA region, the implementation of pilot projects for the development of unused land, such as low-slope 
hilly land, resulted in significant slope-climbing of regional urban construction land (Zhou et al. 2021b; Peng et al. 
2023). 

Used the natural breaks method as basis to classify the BCI values into the four categories of high, middle, low and none 
(Fig. 7e) level, with 9, 33, 86 and 78 cities, respectively. The slope climbing intensity of construction land is mainly 
concentrated in the eastern and southern regions, and only 4 counties in the western region have a medium or high 
level of climbing intensity. Fig. 7f shows that there is only one county with an average ULS value above 20?from 2000 to 
2020, which is Ji County in the HGA. Most of the counties with an average ULS of 10-20? are located in the HGA and 
PHA, and the RMA areas are slightly distributed 

The urban construction land in the MRYR during the research period showed a climbing trend at both the 
geomorphological and county scales, but there were certain differences in certain regions. The slope climbing trend of 
urban construction land in the HGA, PHA, PA, and RMA is the same at two scales, while in the SA, there are certain 
differences in the analysis results at the two scales. On the geomorphological scale, the phenomenon of construction 
land climbing in the SA is not obvious, and construction land tends to expand towards low slope areas; At the county 
level, different towns in the SA exhibit varying degrees of slope climbing at different times. This indicates that regional 
differences at the county level should be fully considered in land use planning and urban construction in the MRYR, 
especially in the SA. 

3.3. Impact of environmental and socio-economic factors on urban construction land expansion 

Based on the SEM results between regional urban construction land area and environmental-economic elements (Tem, 
Pre, PET, and NDVI) and economic factors (Pop and NTL), this study systematically discussed the driving mechanism 
of regional urban construction land expansion at different slopes, slope aspect, and altitude. Regional environmental 
and economic factors (except for Tem) significantly influenced urban construction land area variation at different slope 


gradients (p « 0.05), which presented apparent heterogeneity (Fig. 8). In regions with slopes of 0° to 5°, POP, NTL, 
PET, and NDVI exhibited significant positive direct effects on the variation in regional construction land area. The 


effect of NTL was more significant, whereas that of NDVI was weaker (Fig. 8a). This result indicated that the urban 
expansion process with low slope was controlled by many factors, among which economic development was the 
dominant factor. In regions within the 5?-10? slope range, only NTL significantly affected the urban construction land 


area (p « 0.05) (Fig.8b), which further proved that economic development is an important factor driving the vertical 
expansion of regional construction land. The positive influence of PET on the change in urban 
construction land area was increasingly significant with increasing regional slope, with PET path 


coefficients for 10°-15° and 15°-20° slopes of 0.41 and 0.44, respectively (Fig. 8c-d). For regions within the 20°-25° 
slope range, the effect generated by Pre gradually became more evident (p < 0.05), and the effects of regional Pop and 


NDVI on urban construction land changed from positive to negative (Fig. 8e). Regional climate factors become 
important limiting factor for the expansion of urban construction land to high slope areas. 


In different slope direction gradient regions (Fig. 8f—i), the variation in regional Pop had a significant regulatory effect 
on regional urban construction land area variation, with a considerable effect on the semi-sunny slope region and a 
weaker effect on the sunny slope region. An effect of regional NTL on urban construction land area variation was 
observed in different slope gradient regions, except in the sunny slope region. PET and NDVI had significant direct 
positive effects on the regional urban construction land area on both sunny and semi-sunny slopes (p < 0.05). The 
regulatory effects of PET and NDVI gradually weakened from sunny to semi-sunny slopes. The effects of regional Tem 


and Pre on urban construction land area variation for different slope aspects were weaker. These results indicated 
that the process of construction land expansion on different slope direction gradient were mainly 


regulated by demographic factor. 

The effect of regional Tem on urban construction land area changed from positive to negative as altitude increased. The 
effect of Pre on the urban construction land area at different altitude gradients was not significant. A significant 
regulatory effect of regional Pop on urban construction land area variation only appeared in regions with altitudes less 


than 700 m. As for reiongs with 700-2100 m altitudes, both regional NTL and NDVI displayed significant regulatory 
effects on the regional urban construction land area (p «0.05), which indicated the expansion process of 
construction land in this region were mainly controlled by the regional economic development 
level and greening level. At altitudes below 700 m, PET positively affected the urban construction land area, 
whereas a significant negative effect was observed at the region of altitudes above 700 m. 


4 Discussion 

4.1 Influence of vertical expansion of construction land on regional arable land in the MRYR 

The continuous expansion of urban construction land and the decline in high-quality arable land places great pressure 
on regional food security (Xu et al. 2020). Previous studies have indicated that climbing construction land can 
effectively alleviate farmland encroachment in fragile areas, which is of great significance for the protection of high- 
quality arable land resources in fragile areas (Zhou et al. 2021b). In the present study, the impact of urban construction 
land slope climbing on regional arable land was further investigated by analyzing the area of urban construction land 
encroaching on flat arable land (arable land with a slope < 2° was defined as flat arable land (Zhou et al. 2021b)) 


resources in the MRYR (Fig. 9a). During 2000-2020, the proportion of flat arable land occupied by regional urban 
construction land in high-climbing areas fluctuated significantly, and a larger proportion (50.77%) of flat arable land 
was occupied by regional urban construction in the 2005-2010 period. From 2015 to 2020, the proportion of urban 
construction land encroachment on flat arable land increased, posing a threat to regional food security (grain 
production has declined significantly in recent years). In low-climbing areas, the proportion of urban construction land 
encroachment on flat arable land areas has shown a significant upward trend since 2010, whereas regional grain output 
has shown a significant decline during the same period. Overall, the encroachment of flat arable land resources has 
been alleviated in high-climbing areas but remains serious in low-climbing areas. According to the land use transfer 
matrix of the MRYR during 2000-2020 (Table 3), about 30580.07km? of cropland was converted to other land use 
types, of which 16.48% was transferred to the construction land, and the transferred area from the cropland to the 
construction land was 5040.0km2. As for the construction land, about 5835.20km? was converted from other land use 
types, of which 86.37% came from the cropland. The results of the regional Land use transfer matrix further clarified 
that the expansion of urban construction land poses a threat to regional food security. 

The proportion of urban construction land encroachment on flat arable land in horizontally expanding areas was 
relatively small (3.7396) during 2005-2010. A significant increasing trend was observed in the 2010-2015 period, and a 
downward trend appeared in 2015-2020. Compared with climbing areas, construction land encroachment on fragile 
farmland was not significant in horizontal expansion areas. 

4.2. Impact of construction land climbing on ecological land use 

The occupation of existing ecological land during urban expansion not only reduces the quality of the living 


environment but also easily triggers regional ecological environmental problems. According to previous study 
(Zhou et al. 2021b), we defined the regional forest, grassland and wetland with a slope of x 2? as the flat ecological land 
in this study. The current study revealed the impact of urban construction land expansion on regional 
ecological land by analyzing the area of urban construction land encroaching on flat ecological 


land in the study area during 2000-2020 (Fig. 9). The results showed a significant downward trend in the proportion 
of urban construction land encroaching on flat ecological land in high-climbing areas from 2000 to 2015. The lowest 
proportion of urban construction land occupying fragile ecological land in high-climbing areas was observed in 2015 
(7.4696). The trend in the proportion of urban construction land encroaching on flat ecological land in the low-climbing 
areas was similar to that in the high-climbing areas. 


In the horizontal expansion area, the proportion of flat ecological land occupied by urban construction land was higher 
than that in the other two areas (high- and low-climbing areas) during the study period, exhibiting a significant upward 
trend during 2000-2015 and a slightly decreasing trend during 2015-2020. Overall, the proportion of urban 
construction land encroaching on flat ecological land during slope climbing was lower than that in the horizontal 
expansion area, indicating that urban construction land slope climbing can effectively alleviate the construction land 
encroachment of flat ecological land. 


4.3 Effect of construction land expansion on regional carbon storage in the MRYR 

The spatial and temporal distribution of carbon reserves in the past 20 years in the study area was calculated based on 
the InVEST model. During 2000-2020, the carbon storage in the MRYR increased from 3380.00x106t to 3421.30x10°t, 
and the total carbon storage in the MRYR increased by 27.3x10°t. The spatial heterogeneity of regional carbon reserves 
was also obviously. The lower carbon reserves were concentrated in the SA, the urban area of the PA and the RMA, 
while the higher carbon reserves were mainly concentrated in the forest areas of the central, eastern and southern areas 
in the MRYR. 

The carbon storage per unit area of each land use type in the MRYR was ranked as follows: the Forest (12255 t-km?) > 
the grassland (10161 t-km?) > the cropland (8805 t-km?) > the construction land (6530 t-km?) > the barren (2634 
t-km?) > the water (1163 t-km?). From 2000 to 2010, the carbon stocks of the forest, grassland, water and construction 
land respectively increased by 49.7x106t, 50.5x10°t, 39.1x106t and 18.3x106t, while the carbon stocks of the cropland 
and barren decreased by 84.9x106t and 6.7x106t, respectively. Among various land use types in the MRYR, the 
grassland had the largest carbon storage, accounting for 39.5196 of the total carbon storage, followed by the cropland 
and the forest, accounting for 29.65% and 28.2% of the total carbon storage of the MRYR in 2020. The construction 
land, the barren, and the water region had relatively smaller carbon storage, accounting for only 2.5396, 0.0796, and 
0.05% of the total carbon storage in the MRYR. 

As for different landforms, the carbon storage in the HGA presented relatively highest value (1167.1x10%t in 2020), 
followed by the PHA (exceeding 1070x106t in 2020). During the past 20 years, the carbon storage of the SA, HGA, and 
PHA areas displayed continuous growth trend. However, the carbon storage of the RMA presented fluctuating variation, 
with a decreasing trend from 2000 to 2010 and increasing trend from 2010 to 2020. The carbon storage in the PA 
shown a decreasing trend over past 20 years, and it has decreased to 581.1 x10°t. 

Significant spatial heterogeneity occurs in the influence of urban expansion on regional carbon storage (Fig. 10 g-k). 
The urban expansion process in the PA had the most significant impact on the regional carbon storage. Due to the 
expansion of urban construction land, the carbon storage in the PA decreased by 6.7 x106t during 2000-2020. Similar 
significant influence of the urban expansion on regional carbon storage variation was also observed in the RMA. The 
urban expansion of the SA was relatively weak, which also lead to relatively weak influence on the variation of regional 
carbon storage. As for the county scale (Fig.10), the urban expansion of the metropolitan distribution areas in the PA 
(such as Xi 'an, Taiyuan, etc.) had a more significant impact on the local carbon storage, among which the urban 
expansion of Xi'an led to the decrease of its regional carbon storage by 0.58x106t over past 20 years. There are only 
three counties displayed increasing trend of carbon storage during 2000-2020, namely the Yonghe County 
(0.007x106t), Daning County (0.005x106t), and Baode County (0.002 x106t). 

4.4 Comparison with existing research and outlook for future research 

Urban expansion is characterized by the sustained interaction between vertical and horizontal growth (Zhou et al. 
2021b; Zhang H et al. 2023). However, Existing studies paid much less attention to the slope-climbing phenomenon of 
construction land in the vertical dimension than of the expansion of construction land in the horizontal direction, 
especially in complex landform regions (Zhou et al. 2021b; Zhang H et al. 2023). Therefore, we comprehensively 
explored the vertical slope-climbing characteristics of urban construction land from complex landforms in two spatial 
scales and the fist fine classification-scale study of regional construction land of the MRYR in this study, which thereby 
compensating for the shortcomings of existing research. In addition, previous studies (Li et al. 2022; Zhang H et al. 
2023) pointed out that relatively higher resolution long-time-series LUCC products with (greater than or equal to 30 m 
resolution) can effectively suppress the analytical errors of the evolution of spatiotemporal patterns of vertical urban. In 
this study, we selected a long series of land use data with a resolution of 30 m x 30 m to best reveal the evolutionary 
characteristics of urban construction land with complex landforms, while reducing the analytical error caused by the 
low accuracy of LUCC products. 

In this study, the horizontal expansion and vertical slope climbing of urban construction land in the MRYR coexist and 
displayed obvious spatial heterogeneity, which is similar to the results of existing hotspot studies (such as Guiyang, 
Shen zhen, and the Yangtze River urban agglomeration) and studies in other regions of the Yellow River basin (such as 
the Lanzhou-Xining) (Zhou et al. 2021b; Shi et al. 2022; Peng et al. 2022). It is worth noting that the proportion of 
construction land area located below 5?slopes in most of geomorphic zones ot the MRYR exceeded the average value in 
China, which indicated that the horizontal urban expansion process in the area is in the stage of rapid development. 
Besides, gradual urban construction land expansion to high-slope regions and the slope-climbing phenomenon were 
observed in the areas above slopes of 5? especially in the HGA and RMA region during the study period. Large-scale 
mineral development, GDP growth, population growth and the limited resources of flat land were the potential reasons 
for the vertical expansion of urban construction land in this region (Zhang H et al. 2023). Therefore, scientific and 
effective spatial planning and tailored development models are essential for regional sustainable urban development 
(Zhang H et al. 2023). In addition, some regulatory measures are also necessary for the areas where the construction 
land slope-climbing phenomenon was obviously and occupied more ecological land. 

In this study, we only studied the horizontal expansion and vertical slope-spectrum characteristics, gradient effects and 


driving factors of construction land, and also discussed the Influence of vertical expansion of construction land on 
regional arable land, ecological space and carbon storage. In the future, we plan to conduct long-time-series research 
and more highly detailed studies of the evolutionary characteristics of construction land using higher-spatial-resolution 
and higher-precision data are needed in future. In addition, we should be considered more variables in quantified the 
influence factor of climbing construction land (such as government policies, living customs and historic preservation) 
and make suggestions for the high-quality development of the MRYR, and contribute wisdom to the high-quality 
protection of the Yellow River basin 


5. Conclusions 


Based on the data of urban construction land in different geomorphic areas in the MRYR and the topographic and the 
construction land slope spectrum, the spatial and temporal distribution pattern of urban construction land at different 
scales, especially at the county scale. The structural equation model is used to analyze the influencing factors of regional 
urban construction land area change under different terrain and geomorphic conditions. The study found: 

During the research period, significant spatial heterogeneity was observed in the spatial distribution of natural factors, 
climate, economy, and other factors in the MRYR. A superior natural environment was observed in the southeast and 
plain valleys, which had a concentrated population distribution and improved economic development. The expansion of 
urban construction land in the MRYR has increased significantly since 2000, with the expansion areas mainly 
concentrated in the Fen and Wei River basins and the most significant expansion occurring in the 2005-2010. 
Significant horizontal expansion of urban construction land was observed in the SA and HGA regions. 


The urban construction land in the MRYR has vertically expanded to areas above 5° slopes, especially in 2005-2010. 
The slope spectrum of urban construction land in each geomorphic type area presented a smooth inverted "V" curve, 


with obviously urban construction land slope-climbing appeared in HGA and RMA. 

During 2000-2020, 68.45% of the total number of counties in the MRYR were climbing counties, of which 37.38% were 
high-climbing counties. A higher proportion of climbing counties was located in the RMA region, whereas a lower 
proportion appeared in the HGA region. 

Regional population density and economic development levels strongly influenced urban construction land area 
variability. The influence of the natural environment and socio-economic factors on urban construction land variability 
displayed gradient variation with the regional slope, slope aspect, and altitude. 

The urban expansion process in slope-climbing areas poses a threat to regional food security. However, the regional 
ecological environment of climbing urban construction land improved significantly; therefore, urban construction land 
slope-climbing somewhat alleviates the ecological security pressure exerted by urban expansion. The urban expansion 
of the metropolitan distribution areas in the PA (such as Xi 'an, Taiyuan, etc.) had a more significant impact on the local 
carbon storage. 
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Fig. 1 Location and geographical overview of the study area. (a) Geographical location map of the study area. (b) 
Annual mean precipitation and temperature. (c) Land use area (%) for 2020. 
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Fig. 2 General research framework in this study. 
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Fig. 3 Spatial distribution characteristics of environmental and socio-economic factors for the middle reaches of the 
Yellow River during 2000, 2010, and 2020. (a): Average annual temperature, (b) annual precipitation, (c) annual 
potential evapotranspiration, (d) NDVI, (e) LUCC, (f) population density, and (g) GDP. NDVI, normalized difference 
vegetation index; LUCC, land use and coverage change; GDP, gross domestic product. 
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Fig. 4 Temporal changes in urban construction land in the middle reaches of the Yellow River (MRYR). Spatial 
distribution of urban construction land in 2000 (a), 2010 (b) and 2020 (c). (d) Temporal changes in urban 


construction land. (e) Time change in urban construction land area in each region. (f) Three-dimensional column chart 
of urban construction land area in each region. 
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Fig. 5 Distribution of regional construction land in different geomorphological zones of the middle reaches of the 
Yellow River (MRYR) between 2000 and 2020. (a-b) MRYR; (c-d) sandstorm area (SA); (e-f) loess hilly-gully area 
(HGA); (g-h) loess plateau-hilly area (PHA); (i-j) plain area (PA); (k-l) rocky mountain area (RMA). 
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Fig. 6 Topographic and construction land slope spectrum and the ABCI and ULSC results for construction land in the 
MRYR and each geomorphological zone from 2000 to 2020. Topographic and construction land slope spectrum in 
different regions: (a) MYRY, (b) SA, (c) HGA, (d) PHA, (e) PA, (f) RMA. (g) ABCI and ULSC of the whole MRYR; (h) 
ABCI of each geomorphological zone; (i) ULSC of each geomorphological zone. ABCI, average built-up land climbing 
index; ULSC, upper limited slope angle change; MRYR, middle reaches of the Yellow River; SA, sandstorm area; PHA, 
loess plateau-hilly area; HGA, loess hilly-gully area; PA, plain area; RMA, rocky mountain area. 
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Fig. 7 Spatial distribution of the ABCI, BCI, average ULS and ULSC at the county scale in the middle reaches of the 
Yellow River from 2000 to 2020 (a-d) ABCI. (e) BCI. (f) Average ULS. (g) Cluster distribution of ABCI and ULSC. ABCI, 
average built-up land climbing index; ULSC, upper limited scope angle change. BCI, built-up land climbing index. 
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Fig. 8 Structural equation modeling results for the regional construction land area and environment-economic 


elements. Slope gradients of (a) 0°-5°, (b) 5°-10°, (c) 10°-15°, (d) 15°-20°, and (e) 20°-25°. (f) Shade slope, (g) sunny 
slope, (h) semi-sunny slope, (i) semi-shady slope, (j) altitude gradient below 700 m, (k) altitude gradient of 700-1,400 
m, (1) altitude gradient of 1,400-2,100m, (m) altitude gradient of 2,100-2,800m. 


80r 
- D- high-climbing QS 
- O- low-climbing E "x 
-A horizontal expansion "M ES 
t A b NC "A 
g 60 T g 60 E 
— gee w ae 
204. cdd ae 
g Dp ccr dl 5 A^ 
E inis R Es - O- high-climbing 
g 40 E .n g 4r -© low-climbing 
e ES " e - A horizontal expansior 
Spies u. 
& ü a a 
Ei $ oO... "n 
Š A- i F ~ ~ ^s -9 
420 gr <20 m c g 
Bes. a A E 
~A 


L L L L 1 L 1 L 
2000-2005 2005-2010 2010-2015 2015-2020 2000-2005 2005-2010 2010-2015 2015-2020 


Fig. 9 Area proportion of flat arable and ecological land (slope x2?) occupied by construction land from 2000 to 
2020. (a) Flat arable land. (b) Flat ecological land. 
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Fig.10 Spatial distribution and variation of carbon storage in the middle reaches of the Yellow River (MRYR) from 2000 
to 2020. 


Table 1 Description of the data used in this study. 


Category 


Topography 
data 


Climatic 
data 


Vegetation 
cover 

data 

Land cover 
data 


Economic 
data 


Data type Year 

Digital Elevation / 

Model (DEM) 

Slope / 

Slope aspect / 

Precipitation (Pre) 2000,2010,2020 

Temperature (Tem) 2000-2020 

Potential 

evapotranspiration 2000-2020 

(PET) 

NDVI 2000-2020 

LUCC 2000,2005,2010 
2015,2020 

GDP 2000,2005,2010 
2015,2020 

Population amount SO SU 

(Pop) 

Nighttime light (NTL) 2000-2020 

Grain output 2000-2020 


Resolution 
30mx30m 


30mx30m 
30mx30m 


ikmxikm 


ikmxikm 


ikmxikm 


ikmxikm 


30mx30m 
ikmxikm 
ikmxikm 


ikmxikm 


/ 


Source 

Shuttle Radar Topography Mission 
(SRTM) V3 

calculated by ArcGIS 


calculated by ArcGIS 

Resource and Environment Science and 
Data Center(https://www.resdc.cn/) 
National Tibetan Plateau Data Center 
(https://data.tpdc.ac.cn/) 


National Tibetan Plateau Data Center 
(https://data.tpdc.ac.cn/) 


Geosptial Data Cloud 
(https://www.gscloud.cn/) 


annual China Land Cover Dataset (CLCD) 
(https: //zenodo.org/records/5816591) 
Resource and Environment Science and 
Data Center(https://www.resdc.cn/) 

U.S. Department of Energy's Oak Ridge 
National Laboratory (ORNL) 

LandScan population density dataset 
NOAA (https://ngdc.noaa.gov/) 


National Earth System Science Data 
Center 


Notes: NDVI, normalized difference vegetation index; LUCC, land use and cover change; GDP, gross domestic product 


Table 2 Carbon density of different land use types in the MRYR (/t.hm?) 


Land use type Cabove CBelow 
Cropland 1.02 14.44 
Forset (include shrub) 7.58 20.76 
Grassland 6.71 15.48 
Water (include wetland) O (0) 
Construction land 0.45 o 
Barren 0.23 (0) 


Table 3 Information of regional Land use transfer matrix (km2) 


Cropland 
Forest 
Shrub 
Grassland 
Water 
Barren 


Construction 
land 
Wetland 


Cropland Forest Shrub Grassland 
94975.18 3418.95 6.15 21711.24 
1227.37 62718.94 136.49 185.42 
29.14 769.74 293.87 527.62 
15820.41 8894.35 120.52 105995.56 
140.33 2.15 o 7.91 
145.84 1.02 (0) 3495.84 
49.02 0.08 (0) 1.47 

o o o 0.002 


Water 
383.64 
0.83 
0.02 
76.47 
680.54 
10.30 


180.81 


(0) 


Csoil Cpead 

69.75 2.84 

90.12 4.09 

77.23 2.19 

11.63 (0) 

64.85 (0) 

26.11 (0) 

Barren Construction land Wetland 
20.09 5040.0 0.01 
0.09 29.23 (0) 
0.15 0.15 (0) 
233.09 607.11 0.01 
0.96 101.85 0.01 
618.12 56.86 (0) 
0.63 7079.37 (0) 

o (0) (0) 


Appendix 1 


1. Nighttime light data processing method 


The nighttime lighting data processing and principles are referenced from Chang et al., (2023) and Zheng et al., (2024). The 
details are as follows. 
a. DMSP/OLS data preprocessing 

(D Data extraction 

DMSP/OLS NTL data was extracted from China's administrative boundaries using a mask. 

@ Mutual correction 

The abnormal fluctuation in the pixel DN value of the multisensory images was the main reason for the discontinuity of the 
image data. 

Mutual correction is performed between the sensors to improve image continuity. According to Zheng et al. (2024), A unary 
quadratic equation was selected to establish the correction model, as follows: 

DNc=aDN?+bDN+c 

(1) 

where DN and DNc represent the grey values of the pixels before and after correction, respectively, and a, b, and c are the 
regression parameters. 

(©) Intra-year fusion 

Because of pixel DN value variations at identical locations within multisensor images of the corresponding year and the 
presence of atypical fluctuations, it is imperative to perform image-to-image continuity correction. Using Formula (2), the image 
data obtained from the distinct sensors were corrected, resulting in unique data for each year from 2000 to 2013. 


0, DN% 570, DN, j =0 


(DN;, p +DN}, i) )/2=0, else 

where DN%, DE DNin, j represent the DN values of the ith pixel from sensors x and y, respectively, in the nth year after 
correction according to Formula (1), and DN, i) represents the corrected DN value of the pixel. n = 2000, 2001 ..., 2007. 

(à) Temporal series correction 

Based on China's developmental status, NTL is in an enhanced state. Therefore, assuming a consistently increasing trend in 
the pixel DN values, the DN value of the next year will be higher than that of the previous year, and the DN values of the missing 
or unstable pixels will be replaced with zero. Based on this assumption, the formula for correcting multiyear image data acquired by 
multiple sensors is as follows: 


DN, - (2) 


0, DN a+, —0 
DN;,;-7| DN- DNa+1,>0 N DNq-1,) >DN pi) (3) 
DN,, ;, else 


where DN.;, i, DNq, j, and DNi+z, i represent the DN values of pixel i in the image for years n-/, n, and n+1, respectively. n 
— 2000, 2001 ..., 2013. 

By applying the aforementioned data correction, we obtained a corrected DMSP/OLS NTL dataset for China from 2000 to 
2013. 


b. NPP/VIIRS data preprocessing 

(D Annual image synthesis and cropping 

Monthly NPP/VIIRS NTL data on a global scale between 2012 and 2020 were averaged annually. The annual composite data 
were extracted from the administrative boundaries of China using a mask. 

© Image correction 

By referring to the previous studies (Guan et al., 2021), radiance values lower than 0.3x10? W-cm?-sr'! were regarded as 
noise and were assigned a value of 0 to eliminate the adverse effects of background noise. Apart from this, logarithmic 
transformation was adopted to suppress the sharp changes in pixel values in the urban core area, which was necessary to maintain 
consistency in the numerical ranges between the NPP/VIIRS and DMSP/OLS images (Li et al., 2022). The formula is as follows: 

NPPigg = In(NPP t 1) (4) 

where NPP represents the original radiance of the NPP/VIIRS image, NPPiog is the log-transformed radiance. In practical 
applications, adding constant 1 to avoid negative values. 

By implementing the aforementioned data correction procedures, a rectified NPP/VIIRS NTL dataset encompassing the 
geographic extent spanning the period of China between 2012 and 2020 was obtained. 


c. DMSP/OLS and NPP/VIIRS data consistency correction 
According to Guan et al. (2021) the following formula was used to correct the consistency of DMSP/OLS and NPP/VIIRS 
images for consistency correction. 
p 1-p 
y 7 Ait (s - A) [sic + segs] (9 
where y is the calibration value after consistency correction of NPP/VIIRS and DMSP/OLS. x is the radiation value of 


NPP/VIIRS image after logarithmic transformation. 47, A2, logx01, logx02, hi, h2 and p are the model parameters. According to 
Guan et al. (2021), the model parameters are: -31.5988, 190.6737, 0.2857, 1.9352, -3.113 2, 0.4410, and 0.0290, respectively. 


2.The correction formula of carbon density parameter 
The carbon density will vary depending on climate, soil properties, and land use, and it needs to be corrected. Therefore, the 
formula from A/am et al. (2013) was used as the correction formula for precipitation factor, and the formula from Giardina et al. 
(2000) and Chen et al. (2007) were used as the correction formula for annual average temperature and biomass carbon density. 
Csp-3.3968xP-3996.1 (R?=0.11) 


(6) 
Cap-6.798e990541P (R?=0.70) 

(7) 
Cpr-28xT-398  (R?=0.47,P<0.01) 

(8) 


where Csp is the soil carbon density based on annual precipitation (kg:m?). Csr and Car are the biomass carbon density based 
on annual precipitation and temperature (kg:m?), respectively. P is the average annual precipitation (mm), and T is the average 
annual temperature (°C). 

Soil carbon density and biomass carbon density were calculated for the whole country and the study area, then a correction 
factor was derived according to the following equations. The carbon density for the study area was the product of the national 
carbon density and the correction factor. 


le 
Kgp= Cis 
(9) ^ 
Kgr- NN 
(10) 
Kg-KprXK pp 
(11) o 
-Csr 
a C"sp 
(12) 


where Kap and Kar are the correction coefficients of precipitation factor and temperature factor for biomass carbon density, 
respectively. C'sp and C"sp are the biomass carbon density of the study area and the national scale obtained from annual 
precipitation. C'gr and Csr are the biomass carbon density of the study area and the national Scale based on annual temperature. 
C'sp and C"'sp are the soil carbon density of the study area and the national scale according to the annual temperature. Ks and Ks are 
the Correction coefficient of biomass carbon density and soil carbon density, respectively. 
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